#Analysis of grouse sequence-capture
library(GeneR)
setwd("/Paterson/Datafiles/grouse")

#read in blast m8 table for cDNA against AllContigs
blast.cDNA <- read.table("grouse_analysis/cDNA_vs_scDNA.out",stringsAsFactors=F)
names(blast.cDNA) <- c("Query.id", "Subject.id", "identity", "alignment.length", "mismatches", "gap.openings", "q.start", "q.end", "s.start", "s.end", "e.value", "bit.score")

readFasta("grouse_seq_cap/AllGrouseContigs.fasta",name="contig00065") #loads sequence into buffer

#dir.create("test_files")

blast.cDNA2 <- blast.cDNA[order(blast.cDNA$Query.id),]
#cDNA_101-1 and ssh_477-1 each have >1000 entries

#get each cDNA contig in turn, find sc contig with the top bit scor,
#the next best that doesn't (much) overlap

cDNA.contigs <- unique(blast.cDNA2$Query.id)

tmp.blast <- blast.cDNA2[blast.cDNA2$Query.id %in% cDNA.contigs[42],]
tmp.blast <- tmp.blast[order(tmp.blast$bit.score,decreasing=TRUE),]
#think about adding another column for number of reads

tmp.inc <- rep(FALSE,dim(tmp.blast)[1])
max.q <- max(c(tmp.blast$q.start,tmp.blast$q.end))
tmp.olap <- rep(0,max.q)

for(i in 1:dim(tmp.blast)[1]){
	if(sum(tmp.olap[tmp.blast[i,"q.start"]:tmp.blast[i,"q.end"]])>8) next
	tmp.olap[tmp.blast[i,"q.start"]:tmp.blast[i,"q.end"]] <- 1
	tmp.inc[i] <- TRUE
	}

#export alignments into embl file
system('grep \'>\' grouse_seq_cap/AllGrouseContigs.fasta|sed \'s/>//\'>test_files/sc_contig_info.txt')
sc_contig_info <- readLines("test_files/sc_contig_info.txt")
sc_contig_info2 <- data.frame(sc.contig=sapply(sc_contig_info,function(X){strsplit(X,'[ ]+')[[1]][1]}),
	length=sub('length=','',sapply(sc_contig_info,function(X){strsplit(X,'[ ]+')[[1]][2]})),
	numreads=sub('numreads=','',sapply(sc_contig_info,function(X){strsplit(X,'[ ]+')[[1]][3]})),
	stringsAsFactors=F)

#turn this into function
tmp.blast2 <- tmp.blast[tmp.inc,]
tmp.sc.contigs <- unique(tmp.blast2$Subject.id[order(tmp.blast2$q.start)])
for(i in 1:length(tmp.sc.contigs)){
	tmp.fileout <- paste("test_files/",sub('lcl\\|','',tmp.sc.contigs[i]),".embl",sep="")
	tmp.cont.info <- sc_contig_info2[grep(sub('lcl\\|','',tmp.sc.contigs[i]),sc_contig_info2$sc.contig),][1,]
	tmp.blast.cont <- tmp.blast2[tmp.blast2$Subject.id==tmp.sc.contigs[i],]
	
	writeEmblLine(file=tmp.fileout,code="ID",append=F,text=tmp.cont.info$sc.contig)
	writeEmblLine(file=tmp.fileout,code="CC",
		text=rownames(tmp.cont.info))
	writeEmblLine(file=tmp.fileout,code="FT",header="source", 
		text=paste("1",tmp.cont.info$length,sep=".."),nextfield=F)
	writeEmblLine(file=tmp.fileout,code="FT",
		text="/organism=\"Lagopus lagopus\"",nextfield=F)
	writeEmblLine(file=tmp.fileout,code="FT",
		text="/mol_type=\"genomic DNA\"")
	
	#deal with complement for mRNA markup
	tmp.join <- character(length=dim(tmp.blast.cont)[1])
	for(j in 1:dim(tmp.blast.cont)[1]){
		if(tmp.blast.cont$s.start[j]>tmp.blast.cont$s.end[j]){
			tmp.join[j] <- paste("complement(",tmp.blast.cont$s.start[j],"..",tmp.blast.cont$s.end[j],")",sep="")
			}else{
			tmp.join[j] <- paste(tmp.blast.cont$s.start[j],"..",tmp.blast.cont$s.end[j],sep="")
			}
		}
	if(length(tmp.join)>1){
		writeEmblLine(file=tmp.fileout,code="FT",header="mRNA",
			text= paste("join(", paste(tmp.join,collapse=","),")",sep="") )
		}else{
		writeEmblLine(file=tmp.fileout,code="FT",header="mRNA",
			text= tmp.join )
		}
	readFasta("grouse_seq_cap/AllGrouseContigs.fasta",name=tmp.cont.info$sc.contig) #loads sequence into buffer
	setStrand(0)
	writeEmblSeq(file=tmp.fileout)	
	
	}

#remember >1 exon per sc contig possible when marking CDSs, done

###### CONCATENATE CONTIGS MATCHING EACH EST
#try creating concantenated sequences for GeneWise and ESTWise

tmp.fileout <- paste("test_files/",tmp.blast2$Query.id[1],"_matched.embl",sep="")
writeEmblLine(file=tmp.fileout,code="ID",append=F,text=paste(tmp.blast2$Query.id[1],"_matched",sep=""))

tmp.sc.contigs2 <- sub('lcl\\|','',tmp.sc.contigs)

tmp.pst <- paste(tmp.sc.contigs2,collapse="; ")
tmp.cont.info <- sc_contig_info2[match(tmp.sc.contigs2,sc_contig_info2$sc.contig),]
tmp.cont.info$length <- as.numeric(tmp.cont.info$length)
tmp.cont.info$numreads <- as.numeric(tmp.cont.info$numreads)
writeEmblLine(file=tmp.fileout,code="CC",text=tmp.pst)
writeEmblLine(file=tmp.fileout,code="FT",header="source", 
		text=paste("1",sum(tmp.cont.info$length)+50*(dim(tmp.cont.info)[1]-1),sep=".."),nextfield=F)
writeEmblLine(file=tmp.fileout,code="FT",
		text="/organism=\"Lagopus lagopus\"",nextfield=F)
writeEmblLine(file=tmp.fileout,code="FT",
		text="/mol_type=\"genomic DNA\"",nextfield=F)

#say where the contigs are
tmp.l <- 1
tmp.cont.info$concat.start <- 0
tmp.cont.info$concat.end <- 0
for(i in 1:length(tmp.sc.contigs)){
	tmp.cont.info$concat.start[i] <- tmp.l
	tmp.l <- tmp.l + tmp.cont.info$length[i] +50 -1
	tmp.cont.info$concat.end[i] <- tmp.l
	tmp.l <- tmp.l +1
	writeEmblLine(file=tmp.fileout,code="FT",header="misc_feature",
		text=paste(tmp.cont.info$concat.start[i],"..",tmp.cont.info$concat.end[i]-50,sep=""),nextfield=F)
	#change header to misc_feature?
	
	writeEmblLine(file=tmp.fileout,code="FT",text=paste("/note=\"",tmp.cont.info$sc.contig[i],"\"",sep=""),nextfield=F)
	tmp.blast.cont <- tmp.blast2[tmp.blast2$Subject.id==tmp.sc.contigs[i],]
	tmp.seqi <- strReadFasta("grouse_seq_cap/AllGrouseContigs.fasta",name=tmp.cont.info$sc.contig[i])
	
		if(tmp.blast.cont$s.start[1]>tmp.blast.cont$s.end[1]){ #reverse complement
		tmp.seqi <- strComp(tmp.seqi)
		writeEmblLine(file=tmp.fileout,code="FT",text="/note=\"reverse\"",nextfield=F)
		}else{
		writeEmblLine(file=tmp.fileout,code="FT",text="/note=\"forward\"",nextfield=F)
		}
	
	#add markup for blast hits
	if(tmp.blast.cont$s.start[1]>tmp.blast.cont$s.end[1]) {
		tmp.s.start <- tmp.cont.info$length[i]+1 - tmp.blast.cont$s.end
		tmp.s.end <- tmp.cont.info$length[i] - tmp.blast.cont$s.start +1
		}else{
		tmp.s.end <- tmp.blast.cont$s.end
		tmp.s.start <- tmp.blast.cont$s.start
		}
	tmp.s.end2 <- tmp.cont.info$concat.start[i] + tmp.s.end-1
	tmp.s.start2 <- tmp.cont.info$concat.start[i] + tmp.s.start-1
	
	tmp.tf <- FALSE
	
	for(j in 1:length(tmp.s.start2)){
		if(i == length(tmp.sc.contigs)&j == length(tmp.s.start2)) tmp.tf <- TRUE
		writeEmblLine(file=tmp.fileout,code="FT",header="misc_feature",
			text=paste(tmp.s.start2[j],"..",tmp.s.end2[j],sep=""),nextfield=F)
		#change header to misc_feature?
		
		writeEmblLine(file=tmp.fileout,code="FT",
			text=paste("/note=\"blast hit to query: q.start, ",tmp.blast.cont$q.start[j], 
			"; q.end, ",tmp.blast.cont$q.end[j],
			"; alignment length, ",tmp.blast.cont$alignment.length[j],
			"; mismatches, ",tmp.blast.cont$mismatches[j],
			"; gap openings, ",tmp.blast.cont$gap.openings[j],
			"; bit.score, ",tmp.blast.cont$bit.score[j],"\"",sep=""), nextfield=tmp.tf)
		}
	
	if(i==1){
		tmp.seqout <- tmp.seqi
		}else{
		tmp.seqout <- paste(tmp.seqout,paste(rep("N",50),collapse=""),tmp.seqi,sep="")
		}
	
	
	}
placeString(tmp.seqout,seqno=0)
setStrand(0)
writeEmblSeq(tmp.fileout)
######
#generate fasta files with chicken sequence
grouse_chick1 <- read.table('seq_cap/annotation/grouse_vs_chick.out',stringsAsFactors=F)
names(grouse_chick1) <- c("Query.id", "Subject.id", "identity", "alignment.length", "mismatches", "gap.openings", "q.start", "q.end", "s.start", "s.end", "e.value", "bit.score")
grouse_chick2 <- grouse_chick1[match(unique(grouse_chick1$Query.id),grouse_chick1$Query.id),]

tmp <- readLines('seq_cap/blast_out/grouse600annot.txt')
grouse600annot <- data.frame(grouse.id=rep("xx",length(tmp)),ensembl.id=rep("xx",length(tmp)),stringsAsFactors=F)
grouse600annot$grouse.id <- sapply(X=tmp,FUN=function(X){
	strsplit(X,'\t')[[1]][1]
	})
grouse600annot$ensembl.id <- sapply(X=tmp,FUN=function(X){
	strsplit(X,'\t')[[1]][2]
	})
grouse600annot$uniprot.id <- sapply(X=tmp,FUN=function(X){
	strsplit(X,'\t')[[1]][5]
	})
grouse600annot$best.annot <- factor("none",levels=c("none","ensembl","uniprot"))
grouse600annot$best.annot[nchar(grouse600annot$uniprot.id)==6] <- "uniprot"
grouse600annot$best.annot[grep('ENSGALP',grouse600annot$ensembl.id)] <- "ensembl"

system("grep \'>\' seq_cap/grouse600set.fasta | sed \'s/>//\' > test_files/captured_ids.txt")
cap_ids <- read.table('test_files/captured_ids.txt',stringsAsFactors=FALSE,col.names="cap.id")

library(biomaRt)
chick.mart <- useMart('ensembl',dataset='ggallus_gene_ensembl')
uniprot.mart <- useMart('uniprot_mart',dataset="UNIPROT")
library(GeneR)

#retrieve sequences from biomart, takes ~30min
for(cid in grouse600annot$ensembl.id[grouse600annot$best.annot=="ensembl"]){
	
	tmp.pep.seq <- getSequence(id=cid,type='ensembl_peptide_id',seqType='peptide',mart=chick.mart)
	if(dim(tmp.pep.seq)[1]==0){
		warning(cid," not found\n")
		next
		}
	strWriteFasta(tmp.pep.seq[1,1],file='test_files/chick_peps.fasta',name=tmp.pep.seq[1,2],append=TRUE)
	}

#grouse600annot[grouse600annot$best.annot=="uniprot",]
#then batch retrieval from uniprot.org	

grouse600annot$prot.fa <- "none"
grouse600annot$prot.fa[grouse600annot$best.annot=="uniprot"] <- grouse600annot$uniprot.id[grouse600annot$best.annot=="uniprot"]
grouse600annot$prot.fa[grouse600annot$best.annot=="ensembl"] <- grouse600annot$ensembl.id[grouse600annot$best.annot=="ensembl"]

#generate contigs combined from ESTs
dir.create('/Paterson/Datafiles/grouse/gContigs')

contig.map <- list()
for(gctg in grouse600annot$grouse.id){
	j.out.file <- paste('/Paterson/Datafiles/grouse/gContigs',gctg,"_matched.embl",sep="")
	contig.map[[gctg]] <- join.contigs(contig=gctg,blast.data=blast.cDNA2,contig.info=sc_contig_info2,out.file=j.out.file)
	}
#generate exonerate files
dir.create('/Paterson/Datafiles/grouse/exonerate')

for(gctg in grouse600annot$grouse.id){
	system('rm /Paterson/Datafiles/grouse/test_files/temp.fasta')
	tmp.fasta <- paste('/Paterson/Datafiles/grouse/gContigs/',gctg,"_matched.fasta",sep="")
	tmp.exfile <- paste('/Paterson/Datafiles/grouse/exonerate/',gctg,"_matched.exon",sep="")
	tmp.chickid <- grouse600annot$prot.fa[grouse600annot$grouse.id==gctg]
	if(tmp.chickid=="none") next
	
	tmp.chick.seq <- strReadFasta('/Paterson/Datafiles/grouse/test_files/chick_peps.fasta',name=tmp.chickid)
	strWriteFasta(tmp.chick.seq,file='/Paterson/Datafiles/grouse/test_files/temp.fasta',name=tmp.chickid,append=FALSE)
	
	tmp.cmd <- paste('exonerate  --model protein2genome --showtargetgff TRUE /Paterson/Datafiles/grouse/test_files/temp.fasta',tmp.fasta,'>',tmp.exfile)
	cat(tmp.cmd,'\n')
	system(tmp.cmd)
		}

#add exon annotation
tmp.exon.files <- dir('/Paterson/Datafiles/grouse/exonerate')
dir.create('/Paterson/Datafiles/grouse/embl_exon')

for(exf in 1:length(tmp.exon.files)){#
	tmp.ex.got <- get.exonerate(paste('/Paterson/Datafiles/grouse/exonerate/',tmp.exon.files[exf],sep=""),gene_name=sub('_matched\\.exon','',tmp.exon.files[exf]))
	
	#if no gene prediction in exonerate file
	if(identical(tmp.ex.got,0)) next
	
	tmp.out.file <- paste('/Paterson/Datafiles/grouse/embl_exon/',paste(sub('_matched\\.exon','',tmp.exon.files[exf]),"_exon.embl",sep=""),sep="")
	tmp.in.file <- paste('/Paterson/Datafiles/grouse/gContigs/',paste(sub('_matched\\.exon','',tmp.exon.files[exf]),"_matched.embl",sep=""),sep="")
	
	add.CDS(embl.file=tmp.in.file,cds.list=tmp.ex.got,out.file=tmp.out.file)
	}

gff3.list <- list()
for(exf in 1:length(tmp.exon.files)){#
	tmp.nm <- sub('_matched\\.exon','',tmp.exon.files[exf])
	gff3.list[[tmp.nm]] <- get.exonerate(paste('/Paterson/Datafiles/grouse/exonerate/',tmp.exon.files[exf],sep=""),gene_name=tmp.nm)
	}

#generate SNP data for each matched contig
readNewblerMapper <- function(map.file){
	#note this expects a grepped map file from Newbler
	
	tmp.df <- read.table(map.file,sep="\t",col.names=c('ref.accno','start.pos','end.pos','ref.nucl','var.nucl','depth','freq'),skip=2,stringsAsFactors=FALSE)
	tmp.df$ref.accno <- sub('^>','',tmp.df$ref.accno)
	tmp.df
	}

grouse1hc <- readNewblerMapper('/Paterson/Datafiles/grouse/grouse_seq_cap/grepped/gGrouse1HCDiffs.txt')

#see SNPwrappers files
tst.list <- captureSNPs(contig.map,c("grouse1hc","grouse2hc"),names(contig.map[2]),c("grouse1CO","grouse2CO"),c("/Paterson/Datafiles/grouse/mapping/grouse1.ace","/Paterson/Datafiles/grouse/mapping/grouse2.ace"))
tmp.base <- baseCoverage(tst.list)

#next task, identify whether synonymous or not, and whether in cds or not

#see readAllSnps.R

#check whether any insertions or deletions in exonerate format
exonerate.files <- dir('/Paterson/Datafiles/grouse/exonerate',full.names=TRUE)
tmp.indel <- "#start reading"
for(i in 1:length(exonerate.files)){#
	tmp.indel <- c(tmp.indel,system(paste("grep \'insertions\'",exonerate.files[i]),intern=TRUE))
	}

writeLines(tmp.indel,'grouse_analysis/indel.txt')

tmp.revcomp <- "#start reading"
for(i in 1:length(exonerate.files)){#
	tmp.revcomp <- c(tmp.revcomp,system(paste("grep \'revcomp\'",exonerate.files[i]),intern=TRUE))
	}
#81 exonerate predictions are from the reverse complement strand


#47 frameshifts, need to change these in Artemis
#not sure that masking out 100N gaps is right approach...
# may end up putting everything downstream out of frame
#perhaps change to aaaa...aaaa
#also there are small Ns in sequence
#have to use strTranslate
#probably also have to set up specific masks for frameshifts by
#comparing sequences before and after manual edit with Artemis
#some exonerate predictions are from the reverse complement

#ok done, using coding_snps3.R
# some tidying up to, some repeated entries from newbler mapper
# see hcCoveragelist2[[87]]
#also where 2 different snps found in different populations

#see readAllSnps.R
hcCoveragelist2 <- vector("list",length(contig.map))
names(hcCoveragelist2) <- names(hcCoveragelist)

source('coding_snps3.R')

#for(i.ctg in 1:50){
#	hcCoveragelist2[[i.ctg]] <- try(codingSNPref(hcCoveragelist[[i.ctg]]))
	#names(hcCoveragelist2)[i.ctg] <- names(hcCoveragelist)[i.ctg]
#	}

for(i.ctg in 1:length(hcCoveragelist2)){
	hcCoveragelist2[[i.ctg]] <- try(codingSNPref(hcCoveragelist[[i.ctg]]))
	#names(hcCoveragelist2)[i.ctg] <- names(hcCoveragelist)[i.ctg]
	}


#read in 

grouse_annot <- data.frame(id = names(hcCoveragelist2),stringsAsFactors=FALSE)
grouse_annot$type <- factor(rep("nonimmune",nrow(grouse_annot)),levels=c("nonimmune","immune"))
imm.ids <- read.table('/Paterson/Datafiles/grouse/seq_cap/commands/grouse_imm_ids.txt',stringsAsFactors=FALSE,col.names="id") 
#123 genes, check whether this is the final list

grouse_annot$type[grouse_annot$id %in% imm.ids$id] <- "immune"

grouse_annot$no.cds.snp <- sapply(X=hcCoveragelist2,function(X){
	if(class(X)!="data.frame") return(NA)
	sum(unlist(X$cds))
	})
grouse_annot$no.dn.snp <- sapply(X=hcCoveragelist2,function(X){
	if(class(X)!="data.frame") return(NA)
	sum(unlist(X$dn))
	})
grouse_annot$no.snps <- sapply(X=hcCoveragelist2,function(X){
	if(class(X)!="data.frame") return(NA)
	nrow(X)
	})

#do some QC on SNPs
trimmedHCsnps <- QCsnp(hcCoveragelist2)


##May 2010
#read in AllDiffs file

grouse1All <- readNewblerMapper('/Paterson/Datafiles/grouse/grouse_seq_cap/grepped/gGrouse1AllDiffs.txt')
grouse2All <- readNewblerMapper('/Paterson/Datafiles/grouse/grouse_seq_cap/grepped/gGrouse2AllDiffs.txt')
grouse3All <- readNewblerMapper('/Paterson/Datafiles/grouse/grouse_seq_cap/grepped/gGrouse3AllDiffs.txt')

embl.files <- data.frame(file.name=dir("/Paterson/Datafiles/grouse/embl_exon2",full.name=TRUE),contig = "xx",stringsAsFactors=FALSE)

embl.files$contig <- sapply(embl.files$file.name,function(X){
	sub('_exon.embl','',strsplit(X,"/")[[1]][6])
	})

grouseSNPdfs <- vector('list',length=nrow(embl.files))
names(grouseSNPdfs) <- embl.files$contig
for(dfi in 1:nrow(embl.files)){ #nrow(embl.files)
	grouseSNPdfs[[dfi]] <- try(codingSNPref(
		getAllSNPs(embl.file= embl.files[dfi,1],grouseMapperDataFrames=list(grouse1= grouse1All,grouse2 = grouse2All,grouse3= grouse3All),tsvFiles=c('mapping/grouse1.tsv','mapping/grouse2.tsv','mapping/grouse3.tsv'),tsvLineDataFrames=list(grouse1=grouse1_lines,grouse2=grouse2_lines,grouse3=grouse3_lines))
		,seq.name=embl.files$contig[dfi]))
	if(class(grouseSNPdfs[[dfi]])=="try-error") cat(embl.files$contig[dfi],'failed\n')
	cat(dfi,'  ',embl.files$contig[dfi],'\n')
	}



